"""
Phase 7: Robustness Probes
===========================
1. Investment/GDP control in rate regression (Kopecky-Taylor vs supply-side)
2. Panel cointegration test (spurious regression check)
3. Portfolio equity outlier check
4. Housing with GDP/capita control
"""

import sys
from pathlib import Path

import numpy as np
import pandas as pd
from scipy import stats

PROJECT_DIR = Path("/mnt/c/demographics_capital_flows/asset_returns")
MULTILATERAL_DIR = PROJECT_DIR.parent / "multilateral"
sys.path.insert(0, str(MULTILATERAL_DIR / "src"))
from model import PanelGLS

PROCESSED_DIR = PROJECT_DIR / "data" / "processed"
TABLES_DIR = PROJECT_DIR / "output" / "tables"

OECD_38 = [
    "AUS", "AUT", "BEL", "CAN", "CHL", "COL", "CRI", "CZE", "DNK", "EST",
    "FIN", "FRA", "DEU", "GRC", "HUN", "ISL", "IRL", "ISR", "ITA", "JPN",
    "KOR", "LVA", "LTU", "LUX", "MEX", "NLD", "NZL", "NOR", "POL", "PRT",
    "SVK", "SVN", "ESP", "SWE", "CHE", "TUR", "GBR", "USA",
]

CONTROLS = ['rgdp_growth', 'inflation', 'fiscal_bal_gdp', 'kaopen', 'nfa_gdp_lag']


def stars(p):
    if p < 0.01: return '***'
    if p < 0.05: return '**'
    if p < 0.10: return '*'
    return ''


def run_model(df, dep_var, regressors, label):
    """Run PanelGLS and return results dict."""
    regressors = [r for r in regressors if r in df.columns]
    if dep_var not in df.columns:
        print(f"  [{label}] {dep_var} missing — skipping")
        return None

    sub = df.dropna(subset=[dep_var] + regressors).copy()
    if len(sub) < 50:
        print(f"  [{label}] Insufficient obs ({len(sub)}) — skipping")
        return None

    gls = PanelGLS()
    gls.fit(sub[dep_var].values, sub[regressors].values,
            sub['iso3'].values, sub['year'].values)

    print(f"\n  [{label}]  N={gls.n_obs}, countries={gls.n_countries}, "
          f"R²={gls.r_squared:.4f}, rho={gls.rho:.3f}")

    results = {
        'label': label, 'dep_var': dep_var,
        'n_obs': gls.n_obs, 'n_countries': gls.n_countries,
        'r_squared': gls.r_squared, 'rho': gls.rho,
    }
    for i, name in enumerate(regressors):
        results[f'coef_{name}'] = gls.beta[i]
        results[f'se_{name}'] = gls.se[i]
        results[f'p_{name}'] = gls.pvalues[i]
        sig = stars(gls.pvalues[i])
        print(f"    {name:<30} {gls.beta[i]:>10.4f} ({gls.se[i]:.4f}) {sig}")

    return results


# ══════════════════════════════════════════════════════════════════════
# PROBE 1: Investment/GDP control — Kopecky-Taylor vs supply-side
# ══════════════════════════════════════════════════════════════════════

def probe1_investment_control(df):
    print("\n" + "=" * 70)
    print("PROBE 1: Investment/GDP Control in Safe Rate Regression")
    print("=" * 70)
    print("If Z₁ weakens when controlling for Investment/GDP,")
    print("the rate result is supply-side (investment demand), not KT (bond demand).")
    print("If Z₁ survives, the portfolio channel has independent power.\n")

    demo_vars = ['Z_1', 'Z_2', 'Z_3']
    controls = [c for c in CONTROLS if c in df.columns]

    results = []

    # Check which investment variable is available
    inv_var = None
    for candidate in ['investment_gdp', 'inv_gdp', 'gfcf_gdp', 'ne_gdi_totl_zs']:
        if candidate in df.columns:
            inv_var = candidate
            break

    # If no investment variable, get from multilateral panel
    if inv_var is None:
        print("  No investment/GDP variable in asset_panel. Checking multilateral panel...")
        full = pd.read_csv(MULTILATERAL_DIR / "followup" / "data" / "processed" / "full_panel.csv",
                           usecols=['iso3', 'year', 'gross_investment_gdp'],
                           low_memory=False)
        full = full.dropna(subset=['gross_investment_gdp'])
        # Drop column if it already exists (from prior merge attempt)
        if 'gross_investment_gdp' in df.columns:
            df = df.drop(columns=['gross_investment_gdp'])
        df = df.merge(full, on=['iso3', 'year'], how='left')
        inv_var = 'gross_investment_gdp'
        n_nonnull = df[inv_var].notna().sum()
        n_countries = df.loc[df[inv_var].notna(), 'iso3'].nunique()
        print(f"  Merged {inv_var}: {n_nonnull} non-null, {n_countries} countries")
        if n_nonnull == 0:
            inv_var = None

    if inv_var is None or inv_var not in df.columns or df[inv_var].notna().sum() == 0:
        print("\n  CANNOT RUN PROBE 1: No investment/GDP data available.")
        return results, df

    # Model A: Baseline (no investment control)
    r_base = run_model(df, 'real_bond_10y', demo_vars + controls,
                       "P1-A: Baseline Z → real 10y")
    if r_base:
        results.append(r_base)

    # Model B: With investment/GDP control
    r_inv = run_model(df, 'real_bond_10y', demo_vars + controls + [inv_var],
                      f"P1-B: Z + {inv_var} → real 10y")
    if r_inv:
        results.append(r_inv)

    # Model C: Same for 3m rate
    r_3m_base = run_model(df, 'real_short_3m', demo_vars + controls,
                          "P1-C: Baseline Z → real 3m")
    if r_3m_base:
        results.append(r_3m_base)

    r_3m_inv = run_model(df, 'real_short_3m', demo_vars + controls + [inv_var],
                         f"P1-D: Z + {inv_var} → real 3m")
    if r_3m_inv:
        results.append(r_3m_inv)

    # Summary
    if r_base and r_inv:
        z1_base = r_base['coef_Z_1']
        z1_inv = r_inv['coef_Z_1']
        p_base = r_base['p_Z_1']
        p_inv = r_inv['p_Z_1']
        pct_change = (z1_inv - z1_base) / abs(z1_base) * 100
        print(f"\n  === PROBE 1 VERDICT (10y) ===")
        print(f"  Z₁ without investment:  {z1_base:.3f} (p={p_base:.4f})")
        print(f"  Z₁ with investment:     {z1_inv:.3f} (p={p_inv:.4f})")
        print(f"  Change: {pct_change:+.1f}%")
        if abs(pct_change) > 30 and p_inv > 0.10:
            print(f"  → SUPPLY-SIDE: Z₁ weakens substantially. Rate result is investment-driven.")
        elif abs(pct_change) < 15:
            print(f"  → KOPECKY-TAYLOR: Z₁ robust. Portfolio channel has independent power.")
        else:
            print(f"  → MIXED: Z₁ partially attenuated. Both channels likely operative.")

    return results, df


# ══════════════════════════════════════════════════════════════════════
# PROBE 2: Panel cointegration test
# ══════════════════════════════════════════════════════════════════════

def probe2_cointegration(df):
    print("\n" + "=" * 70)
    print("PROBE 2: Panel Cointegration Test")
    print("=" * 70)
    print("Testing whether Z₁ and real_bond_10y are cointegrated.")
    print("If NOT cointegrated, the level regression may be spurious.\n")

    results = {}

    # Step 1: Panel unit root tests (ADF on each country, then combine)
    # Levin-Lin-Chu style: run ADF per country, combine p-values
    for var_name, var_label in [('Z_1', 'Z₁'), ('real_bond_10y', 'Real 10y yield')]:
        if var_name not in df.columns:
            continue

        adf_stats = []
        countries_tested = 0
        for iso3, grp in df.dropna(subset=[var_name]).groupby('iso3'):
            series = grp.sort_values('year')[var_name].values
            if len(series) < 15:
                continue
            # ADF test: regress Δy on y_{t-1} + Δy_{t-1}
            dy = np.diff(series)
            y_lag = series[:-1]
            if len(dy) < 10:
                continue
            # Simple ADF: regress dy on y_lag
            X = np.column_stack([y_lag[1:], dy[:-1]])
            y = dy[1:]
            if len(y) < 10:
                continue
            try:
                beta = np.linalg.lstsq(X, y, rcond=None)[0]
                resid = y - X @ beta
                se = np.sqrt(np.sum(resid**2) / (len(y) - 2) /
                             np.sum((X[:, 0] - X[:, 0].mean())**2))
                t_stat = beta[0] / se
                adf_stats.append(t_stat)
                countries_tested += 1
            except Exception:
                continue

        if adf_stats:
            mean_t = np.mean(adf_stats)
            # Under H0 (unit root), individual t-stats ~ non-standard dist
            # IPS test: standardized mean approaches N(0,1)
            # Critical value for ADF at 5%: approximately -2.86
            n_reject = sum(1 for t in adf_stats if t < -2.86)
            print(f"  {var_label}: {countries_tested} countries tested")
            print(f"    Mean ADF t-stat: {mean_t:.3f}")
            print(f"    Reject unit root (5%): {n_reject}/{countries_tested} "
                  f"({100*n_reject/countries_tested:.0f}%)")
            if n_reject / countries_tested < 0.3:
                print(f"    → Evidence of UNIT ROOT in {var_label}")
            else:
                print(f"    → {var_label} appears STATIONARY in most panels")
            results[f'unit_root_{var_name}'] = {
                'mean_t': mean_t, 'n_reject': n_reject, 'n_tested': countries_tested
            }

    # Step 2: Engle-Granger cointegration test
    # Run level regression Z₁ → real_bond_10y per country, test residuals for stationarity
    print(f"\n  Engle-Granger cointegration test (Z₁, real_bond_10y):")
    eg_stats = []
    countries_tested = 0
    for iso3, grp in df.dropna(subset=['Z_1', 'real_bond_10y']).groupby('iso3'):
        grp = grp.sort_values('year')
        y = grp['real_bond_10y'].values
        x = grp['Z_1'].values
        if len(y) < 15:
            continue

        # OLS level regression
        X = np.column_stack([np.ones(len(x)), x])
        beta = np.linalg.lstsq(X, y, rcond=None)[0]
        resid = y - X @ beta

        # ADF on residuals
        d_resid = np.diff(resid)
        resid_lag = resid[:-1]
        if len(d_resid) < 10:
            continue
        try:
            X_adf = np.column_stack([resid_lag[1:], d_resid[:-1]])
            y_adf = d_resid[1:]
            b = np.linalg.lstsq(X_adf, y_adf, rcond=None)[0]
            r = y_adf - X_adf @ b
            se = np.sqrt(np.sum(r**2) / (len(y_adf) - 2) /
                         np.sum((X_adf[:, 0] - X_adf[:, 0].mean())**2))
            t_stat = b[0] / se
            eg_stats.append(t_stat)
            countries_tested += 1
        except Exception:
            continue

    if eg_stats:
        mean_eg = np.mean(eg_stats)
        # Engle-Granger critical values are more negative than standard ADF
        # ~-3.34 at 5% for EG with 1 regressor
        n_coint = sum(1 for t in eg_stats if t < -3.34)
        print(f"    Countries tested: {countries_tested}")
        print(f"    Mean EG t-stat: {mean_eg:.3f}")
        print(f"    Reject no-cointegration (5%): {n_coint}/{countries_tested} "
              f"({100*n_coint/countries_tested:.0f}%)")
        if n_coint / countries_tested > 0.5:
            print(f"    → COINTEGRATED: Level regression likely valid")
        elif n_coint / countries_tested > 0.3:
            print(f"    → MIXED: Some evidence of cointegration")
        else:
            print(f"    → NOT COINTEGRATED: Level regression may be SPURIOUS")
        results['cointegration'] = {
            'mean_eg': mean_eg, 'n_coint': n_coint, 'n_tested': countries_tested
        }

    # Step 3: Kao panel cointegration test (pooled residual ADF)
    print(f"\n  Kao-style pooled residual ADF:")
    all_resids = []
    for iso3, grp in df.dropna(subset=['Z_1', 'real_bond_10y']).groupby('iso3'):
        grp = grp.sort_values('year')
        y = grp['real_bond_10y'].values
        x = grp['Z_1'].values
        if len(y) < 10:
            continue
        # Demean within country (FE-style)
        y_dm = y - y.mean()
        x_dm = x - x.mean()
        if np.std(x_dm) < 1e-10:
            continue
        beta = np.sum(x_dm * y_dm) / np.sum(x_dm**2)
        resid = y_dm - beta * x_dm
        all_resids.extend(list(zip(resid[:-1], np.diff(resid))))

    if all_resids:
        resid_lag, d_resid = zip(*all_resids)
        resid_lag = np.array(resid_lag)
        d_resid = np.array(d_resid)
        # Pooled ADF
        b_pool = np.sum(resid_lag * d_resid) / np.sum(resid_lag**2)
        e = d_resid - b_pool * resid_lag
        se_pool = np.sqrt(np.sum(e**2) / (len(e) - 1) / np.sum(resid_lag**2))
        t_pool = b_pool / se_pool
        print(f"    Pooled ADF coefficient: {b_pool:.4f}")
        print(f"    Pooled t-stat: {t_pool:.3f}")
        # Kao critical value at 5%: approximately -1.645 (one-sided)
        if t_pool < -1.645:
            print(f"    → REJECT no-cointegration at 5%. Pooled evidence supports cointegration.")
        else:
            print(f"    → FAIL TO REJECT no-cointegration. Spurious regression concern remains.")
        results['kao_pooled'] = {'b': b_pool, 't': t_pool}

    return results


# ══════════════════════════════════════════════════════════════════════
# PROBE 3: Portfolio equity outlier check
# ══════════════════════════════════════════════════════════════════════

def probe3_portfolio_outliers(df):
    print("\n" + "=" * 70)
    print("PROBE 3: Portfolio Equity/GDP Outlier Check")
    print("=" * 70)

    demo_vars = ['Z_1', 'Z_2', 'Z_3']
    controls = [c for c in CONTROLS if c in df.columns]
    results = []

    if 'port_eq_assets_gdp' not in df.columns:
        print("  port_eq_assets_gdp not in panel — skipping")
        return results

    # Identify outliers
    peq = df.dropna(subset=['port_eq_assets_gdp'])
    print(f"  Full sample: {peq['iso3'].nunique()} countries, "
          f"{len(peq)} obs")

    # Show top countries by mean portfolio equity/GDP
    top = peq.groupby('iso3')['port_eq_assets_gdp'].mean().sort_values(ascending=False)
    print(f"\n  Top 10 countries by mean portfolio equity/GDP:")
    for iso3, val in top.head(10).items():
        n = peq[peq['iso3'] == iso3].shape[0]
        print(f"    {iso3}: {val:.1f}% (N={n})")

    # Model A: Full sample (baseline)
    r_full = run_model(df, 'port_eq_assets_gdp', demo_vars + controls,
                       "P3-A: Full sample Z → port eq/GDP")
    if r_full:
        results.append(r_full)

    # Model B: Drop countries with mean port_eq > 100% of GDP
    high_outliers = top[top > 100].index.tolist()
    print(f"\n  Dropping financial centers (port_eq_assets_gdp > 100% avg): {high_outliers}")
    df_no_fc = df[~df['iso3'].isin(high_outliers)].copy()
    r_no_fc = run_model(df_no_fc, 'port_eq_assets_gdp', demo_vars + controls,
                        "P3-B: Excl financial ctrs → port eq/GDP")
    if r_no_fc:
        results.append(r_no_fc)

    # Model C: Winsorize at p1/p99
    peq_vals = df['port_eq_assets_gdp'].dropna()
    p1, p99 = peq_vals.quantile(0.01), peq_vals.quantile(0.99)
    print(f"\n  Winsorizing at p1={p1:.1f}, p99={p99:.1f}")
    df_wins = df.copy()
    df_wins['port_eq_assets_gdp'] = df_wins['port_eq_assets_gdp'].clip(p1, p99)
    r_wins = run_model(df_wins, 'port_eq_assets_gdp', demo_vars + controls,
                       "P3-C: Winsorized → port eq/GDP")
    if r_wins:
        results.append(r_wins)

    # Model D: Drop top 5 countries
    top5 = top.head(5).index.tolist()
    print(f"\n  Dropping top 5 countries: {top5}")
    df_no5 = df[~df['iso3'].isin(top5)].copy()
    r_no5 = run_model(df_no5, 'port_eq_assets_gdp', demo_vars + controls,
                      "P3-D: Drop top 5 → port eq/GDP")
    if r_no5:
        results.append(r_no5)

    # Summary
    if r_full and r_no_fc:
        print(f"\n  === PROBE 3 VERDICT ===")
        print(f"  Full:        Z₁={r_full['coef_Z_1']:.2f} (p={r_full['p_Z_1']:.4f}), "
              f"R²={r_full['r_squared']:.4f}")
        if r_no_fc:
            print(f"  No fin ctrs: Z₁={r_no_fc['coef_Z_1']:.2f} (p={r_no_fc['p_Z_1']:.4f}), "
                  f"R²={r_no_fc['r_squared']:.4f}")
        if r_wins:
            print(f"  Winsorized:  Z₁={r_wins['coef_Z_1']:.2f} (p={r_wins['p_Z_1']:.4f}), "
                  f"R²={r_wins['r_squared']:.4f}")
        if r_no5:
            print(f"  Drop top 5:  Z₁={r_no5['coef_Z_1']:.2f} (p={r_no5['p_Z_1']:.4f}), "
                  f"R²={r_no5['r_squared']:.4f}")

    return results


# ══════════════════════════════════════════════════════════════════════
# PROBE 4: Housing with GDP per capita control
# ══════════════════════════════════════════════════════════════════════

def probe4_housing_gdppc(df):
    print("\n" + "=" * 70)
    print("PROBE 4: Housing with GDP per Capita Control")
    print("=" * 70)

    results = []
    age_vars = ['old_dep', 'youth_dep']
    housing_controls = ['rgdp_growth', 'real_bond_10y', 'inflation', 'kaopen']
    housing_controls = [c for c in housing_controls if c in df.columns]

    # Check for GDP per capita
    gdppc_var = None
    for candidate in ['gdp_per_capita', 'rgdp_per_capita', 'log_gdp_pc', 'gdppc']:
        if candidate in df.columns:
            gdppc_var = candidate
            break

    if gdppc_var is None:
        print("  No GDP per capita in asset_panel. Checking multilateral panel...")
        full = pd.read_csv(MULTILATERAL_DIR / "followup" / "data" / "processed" / "full_panel.csv",
                           usecols=['iso3', 'year', 'gdp_pc_ppp'],
                           low_memory=False)
        full = full.dropna(subset=['gdp_pc_ppp'])
        full['log_gdppc'] = np.log(full['gdp_pc_ppp'].clip(lower=1))
        if 'log_gdppc' in df.columns:
            df = df.drop(columns=['log_gdppc'])
        df = df.merge(full[['iso3', 'year', 'log_gdppc']], on=['iso3', 'year'], how='left')
        gdppc_var = 'log_gdppc'
        n_nonnull = df[gdppc_var].notna().sum()
        n_countries = df.loc[df[gdppc_var].notna(), 'iso3'].nunique()
        print(f"  Merged log_gdppc: {n_nonnull} non-null, {n_countries} countries")
        if n_nonnull == 0:
            gdppc_var = None

    if gdppc_var is None or gdppc_var not in df.columns or df[gdppc_var].notna().sum() == 0:
        print("  CANNOT RUN PROBE 4: No GDP per capita data.")
        return results, df

    # Model A: Baseline age decomposition (no GDP/capita)
    r_base = run_model(df, 'd_rhpi', age_vars + housing_controls,
                       "P4-A: age ratios → Δreal HPI")
    if r_base:
        results.append(r_base)

    # Model B: With GDP per capita
    r_gdppc = run_model(df, 'd_rhpi', age_vars + housing_controls + [gdppc_var],
                        f"P4-B: age ratios + {gdppc_var} → Δreal HPI")
    if r_gdppc:
        results.append(r_gdppc)

    # Model C: With lagged old_dep (5 year)
    df_lag = df.copy()
    df_lag['old_dep_lag5'] = df_lag.groupby('iso3')['old_dep'].shift(5)
    df_lag['youth_dep_lag5'] = df_lag.groupby('iso3')['youth_dep'].shift(5)
    lag_age = ['old_dep_lag5', 'youth_dep_lag5']
    r_lag = run_model(df_lag, 'd_rhpi', lag_age + housing_controls,
                      "P4-C: lagged age ratios → Δreal HPI")
    if r_lag:
        results.append(r_lag)

    # Model D: With both GDP/capita and lagged age
    r_both = run_model(df_lag, 'd_rhpi', lag_age + housing_controls + [gdppc_var],
                       f"P4-D: lagged age + {gdppc_var} → Δreal HPI")
    if r_both:
        results.append(r_both)

    # Summary
    print(f"\n  === PROBE 4 VERDICT ===")
    for r in results:
        old_key = [k for k in ['old_dep', 'old_dep_lag5'] if f'coef_{k}' in r]
        yth_key = [k for k in ['youth_dep', 'youth_dep_lag5'] if f'coef_{k}' in r]
        old_str = f"{r[f'coef_{old_key[0]}']:.1f} (p={r[f'p_{old_key[0]}']:.3f})" if old_key else "—"
        yth_str = f"{r[f'coef_{yth_key[0]}']:.1f} (p={r[f'p_{yth_key[0]}']:.3f})" if yth_key else "—"
        print(f"  {r['label']:<45} old_dep={old_str}, youth_dep={yth_str}")

    return results, df


# ══════════════════════════════════════════════════════════════════════
# Main
# ══════════════════════════════════════════════════════════════════════

def main():
    print("=" * 70)
    print("PHASE 7: Robustness Probes")
    print("=" * 70)

    df = pd.read_csv(PROCESSED_DIR / "asset_panel.csv")
    print(f"Panel: {df['iso3'].nunique()} countries, {len(df)} obs\n")

    all_md = ["# Phase 7: Robustness Probes\n"]

    # Probe 1
    p1_results, df = probe1_investment_control(df)

    # Probe 2
    p2_results = probe2_cointegration(df)

    # Probe 3
    p3_results = probe3_portfolio_outliers(df)

    # Probe 4
    p4_results, df = probe4_housing_gdppc(df)

    # Build summary table
    print("\n\n" + "=" * 70)
    print("BUILDING SUMMARY TABLE")
    print("=" * 70)

    all_md.append("## Probe 1: Investment/GDP Control (Kopecky-Taylor vs Supply-Side)\n")
    if p1_results:
        all_md.append("| Model | Z₁ | p | R² | N |")
        all_md.append("|---|---|---|---|---|")
        for r in p1_results:
            z1 = r.get('coef_Z_1', float('nan'))
            p = r.get('p_Z_1', float('nan'))
            if not np.isnan(z1):
                all_md.append(f"| {r['label']} | {z1:.3f}{stars(p)} | {p:.4f} "
                              f"| {r['r_squared']:.3f} | {r['n_obs']} |")
    else:
        all_md.append("No investment/GDP data available.\n")

    all_md.append("\n## Probe 2: Panel Cointegration\n")
    if 'cointegration' in p2_results:
        coint = p2_results['cointegration']
        all_md.append(f"- Engle-Granger: {coint['n_coint']}/{coint['n_tested']} countries "
                      f"reject no-cointegration at 5%")
        all_md.append(f"- Mean EG t-stat: {coint['mean_eg']:.3f}")
    if 'kao_pooled' in p2_results:
        kao = p2_results['kao_pooled']
        all_md.append(f"- Kao pooled ADF: t={kao['t']:.3f} "
                      f"({'reject' if kao['t'] < -1.645 else 'fail to reject'} at 5%)")
    for var in ['Z_1', 'real_bond_10y']:
        key = f'unit_root_{var}'
        if key in p2_results:
            ur = p2_results[key]
            all_md.append(f"- Unit root ({var}): {ur['n_reject']}/{ur['n_tested']} reject "
                          f"({100*ur['n_reject']/ur['n_tested']:.0f}%)")

    all_md.append("\n## Probe 3: Portfolio Equity Outlier Check\n")
    if p3_results:
        all_md.append("| Model | Z₁ | p | R² | N | Countries |")
        all_md.append("|---|---|---|---|---|---|")
        for r in p3_results:
            z1 = r.get('coef_Z_1', float('nan'))
            p = r.get('p_Z_1', float('nan'))
            if not np.isnan(z1):
                all_md.append(f"| {r['label']} | {z1:.2f}{stars(p)} | {p:.4f} "
                              f"| {r['r_squared']:.4f} | {r['n_obs']} | {r['n_countries']} |")

    all_md.append("\n## Probe 4: Housing with GDP per Capita\n")
    if p4_results:
        all_md.append("| Model | old_dep | p | youth_dep | p | R² | N |")
        all_md.append("|---|---|---|---|---|---|---|")
        for r in p4_results:
            for ok in ['old_dep', 'old_dep_lag5']:
                if f'coef_{ok}' in r:
                    old_c, old_p = r[f'coef_{ok}'], r[f'p_{ok}']
                    break
            else:
                old_c, old_p = float('nan'), float('nan')
            for yk in ['youth_dep', 'youth_dep_lag5']:
                if f'coef_{yk}' in r:
                    yth_c, yth_p = r[f'coef_{yk}'], r[f'p_{yk}']
                    break
            else:
                yth_c, yth_p = float('nan'), float('nan')
            if not np.isnan(old_c):
                all_md.append(f"| {r['label']} | {old_c:.1f}{stars(old_p)} | {old_p:.3f} "
                              f"| {yth_c:.1f}{stars(yth_p)} | {yth_p:.3f} "
                              f"| {r['r_squared']:.3f} | {r['n_obs']} |")

    out_path = TABLES_DIR / "robustness_probes.md"
    out_path.write_text('\n'.join(all_md))
    print(f"\nSaved: {out_path}")

    print("\n" + "=" * 70)
    print("Phase 7 complete.")
    print("=" * 70)


if __name__ == "__main__":
    main()
